home *** CD-ROM | disk | FTP | other *** search
/ PsL Monthly 1993 December / PSL Monthly Shareware CD-ROM (December 1993).iso / prgmming / dos / c / prime2.exe / PRIME2.C < prev    next >
C/C++ Source or Header  |  1992-09-06  |  11KB  |  362 lines
  1. /******************************************************************************
  2.  * prime2.c
  3.  *
  4.  * (C) 1992 L.I.Kirby    CIS: 70734,126
  5.  *
  6.  * 6th September 1992
  7.  *
  8.  * A couple of algorithms to list and count primes up to reasonably large
  9.  * totals (say 2^32).
  10.  *
  11.  * Both algorithms make use of rolling totals to generate multiples of the base
  12.  * primes, i.e. the primes < sqrt(targetprime)
  13.  *****************************************************************************/
  14.  
  15.  
  16. #include <stdio.h>
  17. #include <ctype.h>
  18.  
  19.  
  20. #define ROLLSIZE    (3405)      /* Big enough to hold all of the base primes */
  21. #define SIEVESIZE   (16000)     /* Must be > largest base prime/2 */
  22. #define SENTINEL    (~0L)       /* Largest possible value for heap endstop */
  23.  
  24. typedef unsigned long LARGE_T;  /* Variables up to targetprime */
  25. typedef unsigned int  SMALL_T;  /* Variables less than sqrt(targetprime) */
  26.  
  27. typedef int bool;
  28. #define FALSE   (0)
  29. #define TRUE    (!0)
  30.  
  31.  
  32.         /* Function prototypes */
  33.  
  34. int     main(int, char **);
  35. static long prime_heap(LARGE_T, bool);
  36. static void downheap(void);
  37. static long prime_sieve(LARGE_T, bool);
  38.  
  39.  
  40.         /* Statics */
  41.  
  42. static LARGE_T S_roll[ROLLSIZE+1];
  43. static SMALL_T S_increm[ROLLSIZE+1];
  44.  
  45.  
  46. /**********************************  main()  **********************************
  47.  * Decode command options and call prime routine.
  48.  */
  49.  
  50. int main(argc, argv)
  51. int     argc;
  52. char    **argv;
  53.  
  54. {
  55.     int     i;
  56.     char    *sptr;
  57.     enum { AL_HEAP, AL_SIEVE } algorithm;
  58.     LARGE_T targetprime;
  59.     bool    printflag;
  60.     long    primecount;
  61.  
  62.     algorithm = AL_SIEVE;
  63.     targetprime = 0L;
  64.     printflag = FALSE;
  65.  
  66.     for (i = 1; i < argc; i++) {
  67.          sptr = argv[i];
  68.         if (*sptr == '-' || *sptr == '/') {
  69.             while (*++sptr) switch (toupper(*sptr)) {
  70.             case 'P':
  71.                 printflag = TRUE;
  72.                 break;
  73.             case 'N':
  74.                 printflag = FALSE;
  75.                 break;
  76.             case 'S':
  77.                 algorithm = AL_SIEVE;
  78.                 break;
  79.             case 'H':
  80.                 algorithm = AL_HEAP;
  81.                 break;
  82.             default:
  83.                 break;
  84.             }
  85.         } else
  86.             (void) sscanf(sptr, "%lu", &targetprime);
  87.     }
  88.  
  89.     if (!targetprime) {
  90.         puts("\nUsage: prime2 [-options] target");
  91.         puts("target:   Search up to but excluding target");
  92.         puts("Options:  P - print primes");
  93.         puts("          N - do not print primes (default)");
  94.         puts("          H - use heap method");
  95.         puts("          S - use sieve method (default)");
  96.         return(1);
  97.     }
  98.  
  99.     switch (algorithm) {
  100.     case AL_HEAP:
  101.         primecount = prime_heap(targetprime, printflag);
  102.         break;
  103.     case AL_SIEVE:
  104.         primecount = prime_sieve(targetprime, printflag);
  105.         break;
  106.     }
  107.  
  108.     printf("\n\nNumber of primes found is %ld\n", primecount);
  109.  
  110.     return(0);
  111. }
  112.  
  113.  
  114. /*******************************  prime_heap()  *******************************
  115.  * This method maintains a running multiple (the roll) for each base prime.
  116.  * These are maintained in a heap structure so that the smallest roll value
  117.  * i.e. next non-prime is easily determined. Even numbers are skipped.
  118.  *
  119.  * Although slower than the sieve method, this method uses less memory and
  120.  * still avoids the dreaded divides.
  121.  */
  122.  
  123.  
  124. static long prime_heap(targetprime, printflag)
  125. LARGE_T targetprime;    /* Maximum search number */
  126. bool    printflag;      /* TRUE is primes should be printed */
  127.  
  128.  
  129. {
  130.     register LARGE_T num; /* The current number being tested for primeness */
  131.     long    nextnonprime; /* The next smallest number which has a factor */
  132.     long    primecount;   /* The number of primes found so far */
  133.      bool    gtsqrt;       /* Set to TRUE once we go beyond sqrt*targetprime) */
  134.     SMALL_T rollhwm;      /* High watermark for roll arrays */
  135.  
  136.  
  137.     printf("\nPrime Heap prime search up to %lu\n", targetprime);
  138.  
  139. /* Initialise array with sentinels for heap */
  140.  
  141.     for (num = 1; num <= ROLLSIZE; num++)
  142.         S_roll[num] = SENTINEL;
  143.  
  144. /* Deal with the special case of 2 */
  145.  
  146.     primecount = 0;
  147.     if (targetprime > 2) {
  148.         primecount++;
  149.         if (printflag)
  150.             printf("\n       2");
  151.     }
  152.  
  153. /* The main routine */
  154.  
  155.     rollhwm = 1;        /* Starting at 1 simplifies the heap routines */
  156.     num = 3;
  157.     nextnonprime = 9;
  158.     gtsqrt = FALSE;
  159.     for (;;) { /* Here all odd numbers between num and nextnonprime are prime */
  160.         do {
  161.             if (num >= targetprime)     /* FINISHED! */
  162.                 return(primecount);
  163.  
  164.             primecount++;               /* Found a prime, count and print it */
  165.             if (printflag)
  166.                 printf(" %7lu", num);
  167.  
  168.             /* If num < sqrt(targetprime) then add to end of heap */
  169.  
  170.             if (!gtsqrt) {
  171.                 if (num * num >= targetprime)
  172.                     gtsqrt = TRUE;
  173.                 else {
  174.                     if (rollhwm >= ROLLSIZE) {
  175.                         printf("\nROLLSIZE too small. Maximum target is %lu\n",
  176.                                                                      num * num);
  177.                         exit(1);
  178.                     }
  179.                     S_roll[rollhwm] = num * num;/* Smallest we need look at -
  180.                                                    guaranteed to be greater than
  181.                                                    any current heap entry */
  182.                     S_increm[rollhwm] = 2 * num;/* Skip over even values */
  183.                     rollhwm++;
  184.                 }
  185.             }
  186.  
  187.             num += 2;
  188.         } while (num < nextnonprime);
  189.  
  190.         /* We've reached nextnonprime. Search for next prime */
  191.  
  192.         do {
  193.             do {
  194.                 S_roll[1] += S_increm[1];   /* Roll to next odd multiple */
  195.                 downheap();                 /* Repair heap */
  196.             } while (num == S_roll[1]);
  197.         } while ((num += 2) == S_roll[1]);
  198.  
  199.          nextnonprime = S_roll[1];
  200.     }
  201. }
  202.  
  203.  
  204. /********************************  downheap()  ********************************
  205.  * Maintains a heap such that:
  206.  * S_roll[n] <= S_roll[2*n] and S_roll[n] < S_roll[2*n+1]
  207.  * i.e. S_roll[1] is the smallest element in the heap.
  208.  *
  209.  * This routine repairs a heap where only the first element is out of order.
  210.  * It expectc unused array elements to contain SENTINEL
  211.  */
  212.  
  213. static void downheap()
  214.  
  215. {
  216.     register int i, j;
  217.     SMALL_T increm;
  218.     LARGE_T roll;
  219.  
  220.     i = 1;
  221.     increm = S_increm[1];
  222.     roll = S_roll[1];
  223.     while ((j = 2 * i) < ROLLSIZE) {
  224.         if (S_roll[j] > S_roll[j+1])
  225.             j++;
  226.  
  227.         if (roll <= S_roll[j])
  228.             break;
  229.  
  230.         S_roll[i] = S_roll[j];
  231.         S_increm[i] = S_increm[j];
  232.  
  233.         i = j;
  234.     }
  235.     S_increm[i] = increm;
  236.     S_roll[i] = roll;
  237. }
  238.  
  239.  
  240. /******************************  prime_sieve()  *******************************
  241.  * This method maintains a running multiple (the roll) for each base prime
  242.  * similar to prime_heap(). However this method uses the multiples to help
  243.  * generate a sieve in sections small enough to fit in memory.
  244.  *
  245.  * S_sieve[] must be at least sqrt(targetprime) in size but larger values
  246.  * will generally be more efficient. Even numbers are ignored.
  247.  */
  248.  
  249.  
  250. static long prime_sieve(targetprime, printflag)
  251. LARGE_T targetprime;    /* Maximum search number */
  252. bool    printflag;      /* TRUE is primes should be printed */
  253.  
  254.  
  255. {
  256.     register SMALL_T rp, sp;
  257.     LARGE_T roll, prime;
  258.     SMALL_T sectlen;
  259.     LARGE_T sectstart, sectend;
  260.     SMALL_T rollhwm;      /* High watermark for roll arrays */
  261.     long    primecount;   /* The number of primes found so far */
  262.     bool    gtsqrt;       /* Set to TRUE once we go beyond sqrt*targetprime) */
  263.     static char S_sieve[SIEVESIZE];
  264.  
  265.      printf("\nSieve prime search up to %lu\n", targetprime);
  266.  
  267.     if (targetprime > (LARGE_T) 4 * SIEVESIZE * SIEVESIZE) {
  268.         printf("\nSIEVESIZE too small. Maximum target is about %lu\n",
  269.